home *** CD-ROM | disk | FTP | other *** search
/ Technotools / Technotools (Chestnut CD-ROM)(1993).ISO / lang_oth / ludecomp / linsys.for
Text File  |  1989-05-14  |  8KB  |  278 lines
  1.       SUBROUTINE LINSYS(A,B,X,N,IORDER)
  2. C
  3. C   SOLVES THE LINEAR SYSTEM AX=B USING LU DECOMPOSITION
  4. C   THIS ROUTINE ASSUMES THAT A IS STORED IN THE MAIN PROGRAM AS A
  5. C   ONE-DIMENSIONAL ARRAY OF DIMENSION N*N, SO THE TRANSPOSE OF A
  6. C   MUST BE TAKEN BEFORE CALLING LUDCMP AND SOLVLU.
  7. C   RETURNS:
  8. C   X       THE SOLUTION TO THE LINEAR SYSTEM
  9. C   A       REPLACED BY ITS LU DECOMPOSITION IN COMPACT FORM
  10. C   B       REORDERED TO CORRESPOND TO ROW INTERCHANGES PERFORMED ON A
  11. C   IORDER  THE ORDER OF THE ELEMENTS OF B
  12. C
  13.       DIMENSION A(N,N),B(N),X(N)
  14.       DIMENSION IORDER(N)
  15. C
  16. C   FIND TRANSPOSE OF A SO THAT THE 2 DIMENSIONAL REFERENCES TO A
  17. C   CORRESPOND TO THE 1 DIMENSIONAL DECLARATION IN THE MAIN PROGRAM
  18. C
  19.       DO 10 I=1,N-1
  20.          DO 20 J=I+1,N
  21.             A(I,J)=A(J,I)
  22.    20    CONTINUE
  23.    10 CONTINUE
  24. C
  25.       CALL LUDCMP(A,N,N,IORDER)
  26.       CALL SOLVLU(A,B,X,N,N,IORDER)
  27.       RETURN
  28.       END
  29. C
  30.       SUBROUTINE LUDCMP(A,N,NDIM,IORDER)
  31. C
  32. C   FROM "APPLIED NUMERICAL ANALYSIS" BY GERALD AND WHEATLEY,
  33. C   THIRD EDITION (ADDISON-WESLEY, 1984).
  34. C
  35. C ---------------------------------------------------------------
  36. C
  37. C   THIS SUBROUTINE COMPUTES THE L AND U TRIANGULAR MATRICES
  38. C   EQUIVALENT TO THE A MATRIX, SUCH THAT LU=A.  THESE MATRICES
  39. C   ARE RETURNED IN THE SPACE OF A, IN COMPACT FORM.  THE U
  40. C   MATRIX HAS ONES ON ITS DIAGONAL.  PARTIAL PIVOTING IS USED TO
  41. C   GIVE MAXIMUM VALUED ELEMENTS ON THE DIAGONAL OF L.  THE ORDER
  42. C   OF THE ROWS AFTER PIVOTING IS RETURNED IN THE INTEGER VECTOR
  43. C   IORDER.  THIS SHOULD BE USED TO REORDER R.H.S. VECTORS BEFORE
  44. C   SOLVING THE SYSTEM AX=B.
  45. C
  46. C ---------------------------------------------------------------
  47. C
  48. C   PARAMETERS:
  49. C
  50. C   A         THE N X N MATRIX OF COEFFICIENTS
  51. C   N         THE NUMBER OF EQUATIONS
  52. C   NDIM      THE FIRST DIMENSION OF A IN THE CALLING PROGRAM
  53. C   IORDER    INTEGER VECTOR HOLDING ROW ORDER AFTER PIVOTING
  54. C
  55. C   THIS ROUTINE CALLS A SUBROUTINE APVT TO LOCATE THE PIVOT ROW AND
  56. C   MAKE INTERCHANGES.
  57. C
  58. C ---------------------------------------------------------------
  59. C
  60.       DIMENSION A(NDIM,N)
  61.       DIMENSION IORDER(N)
  62. C
  63. C ---------------------------------------------------------------
  64. C
  65. C   ESTABLISH INITIAL ORDERING IN ORDER VECTOR
  66. C
  67.       DO 10 I=1,N
  68.          IORDER(I)=I
  69.    10 CONTINUE
  70. C
  71. C   DO PIVOTING FOR FIRST COLUMN BY CALL TO SUBROUTINE APVT
  72. C
  73.       CALL APVT(A,N,NDIM,IORDER,1)
  74. C
  75. C ---------------------------------------------------------------
  76. C
  77. C   IF PIVOT ELEMENT VERY SMALL, PRINT ERROR MESSAGE AND RETURN
  78. C
  79.       IF (ABS(A(1,1)).LT.1.0E-5) THEN
  80.          PRINT 100
  81.          RETURN
  82.       END IF
  83. C
  84. C   NOW COMPUTE ELEMENTS FOR FIRST ROW OF U
  85. C
  86.       DO 20 KCOL=2,N
  87.          A(1,KCOL)=A(1,KCOL)/A(1,1)
  88.    20 CONTINUE
  89. C
  90. C ---------------------------------------------------------------
  91. C
  92. C   COMPLETE THE COMPUTING OF L AND U ELEMENTS.  THE GENERAL PLAN IS TO
  93. C   COMPUTE A COLUMN OF L'S, THEN CALL APVT TO INTERCHANGE ROWS, AND THEN
  94. C   GET A ROW OF U'S.
  95. C
  96.       NM1=N-1
  97.       DO 80 JCOL=2,NM1
  98. C
  99. C   FIRST COMPUTE A COLUMN OF L'S
  100. C
  101.       JM1=JCOL-1
  102.       DO 50 IROW=JCOL,N
  103.          SUM=0
  104.          DO 40 KCOL=1,JM1
  105.             SUM=SUM+A(IROW,KCOL)*A(KCOL,JCOL)
  106.    40    CONTINUE
  107.          A(IROW,JCOL)=A(IROW,JCOL)-SUM
  108.    50 CONTINUE
  109. C
  110. C ---------------------------------------------------------------
  111. C
  112. C   NOW INTERCHANGE ROWS IF NEED BE, THEN TEST FOR TOO SMALL PIVOT
  113. C
  114.       CALL APVT(A,N,NDIM,IORDER,JCOL)
  115.       IF (ABS(A(JCOL,JCOL)).LT.1.0E-5) THEN
  116.          PRINT 100
  117.          RETURN
  118.       END IF
  119. C
  120. C   NOW WE GET A ROW OF U'S
  121. C
  122.       JP1=JCOL+1
  123.       DO 70 KCOL=JP1,N
  124.          SUM=0
  125.          DO 60 IROW=1,JM1
  126.             SUM=SUM+A(JCOL,IROW)*A(IROW,KCOL)
  127.    60    CONTINUE
  128.          A(JCOL,KCOL)=(A(JCOL,KCOL)-SUM)/A(JCOL,JCOL)
  129.    70 CONTINUE         
  130.    80 CONTINUE
  131. C
  132. C ---------------------------------------------------------------
  133. C
  134. C   STILL NEED TO GET LAST ELEMENT IN L MATRIX
  135. C
  136.       SUM=0
  137.       DO 90 KCOL=1,NM1
  138.          SUM=SUM+A(N,KCOL)*A(KCOL,N)
  139.    90 CONTINUE
  140.       A(N,N)=A(N,N)-SUM
  141.       RETURN
  142. C
  143.   100 FORMAT(/,' VERY SMALL PIVOT ELEMENT INDICATES A NEARLY',
  144.      +       ' SINGULAR MATRIX.')
  145.       END
  146. C
  147.       SUBROUTINE APVT(A,N,NDIM,IORDER,JCOL)
  148. C
  149. C   FROM "APPLIED NUMERICAL ANALYSIS" BY GERALD AND WHEATLEY,
  150. C   THIRD EDITION (ADDISON-WESLEY, 1984).
  151. C
  152. C ---------------------------------------------------------------
  153. C
  154. C   THIS SUBROUTINE FINDS THE LARGEST ELEMENT FOR PIVOT IN JCOL OF
  155. C   MATRIX A, PERFORMS INTERCHANGES OF ELEMENTS IN A AND ALSO
  156. C   INTERCHANGES THE ELEMENTS IN THE ORDER VECTOR.
  157. C
  158. C ---------------------------------------------------------------
  159. C
  160. C   PARAMETERS:
  161. C
  162. C   A         MATRIX OF COEFFICIENTS WHOSE ROWS ARE TO BE INTERCHANGED
  163. C   N         NUMBER OF EQUATIONS
  164. C   NDIM      FIRST DIMENSION OF A IN THE MAIN PROGRAM
  165. C   IORDER    INTEGER VECTOR TO HOLD ROW ORDERING
  166. C
  167. C ---------------------------------------------------------------
  168. C
  169.       DIMENSION A(NDIM,N)
  170.       DIMENSION IORDER(N)
  171. C
  172. C ---------------------------------------------------------------
  173. C
  174. C  FIND PIVOT ROW, CONSIDERING ONLY THE ELEMENTS ON AND BELOW DIAGONAL
  175. C
  176.       IPVT=JCOL
  177.       BIG=ABS(A(JCOL,JCOL))
  178.       JP1=JCOL+1
  179.       DO 10 IROW=JP1,N
  180.          ANEXT=ABS(A(IROW,JCOL))
  181.          IF (ANEXT.GT.BIG) IPVT=IROW
  182.    10 CONTINUE
  183. C
  184. C ---------------------------------------------------------------
  185. C
  186. C   NOW INTERCHANGE ROW ELEMENTS IN THE ROW WHOSE NUMBER EQUALS JCOL WITH
  187. C   THE PIVOT ROW UNLESS PIVOT ROW IS JCOL.
  188. C
  189.       IF (IPVT.EQ.JCOL) THEN
  190.          RETURN
  191.       END IF
  192.       DO 20 KCOL=1,N
  193.          SAVE=A(JCOL,KCOL)
  194.          A(JCOL,KCOL)=A(IPVT,KCOL)
  195.          A(IPVT,KCOL)=SAVE
  196.    20 CONTINUE
  197. C
  198. C ---------------------------------------------------------------
  199. C
  200. C   NOW SWITCH ELEMENTS IN THE ORDER VECTOR
  201. C
  202.       ISAVE=IORDER(JCOL)
  203.       IORDER(JCOL)=IORDER(IPVT)
  204.       IORDER(IPVT)=ISAVE
  205.       RETURN
  206.       END
  207. C
  208.       SUBROUTINE SOLVLU(RLU,B,X,N,NDIM,IORDER)
  209. C
  210. C   FROM "APPLIED NUMERICAL ANALYSIS" BY GERALD AND WHEATLEY,
  211. C   THIRD EDITION (ADDISON-WESLEY, 1984).
  212. C
  213. C ---------------------------------------------------------------
  214. C
  215. C   THIS SUBROUTINE IS USED TO FIND THE SOLUTION TO A SYSTEM OF EQUATIONS,
  216. C   AX=B, AFTER THE LU EQUIVALENT OF A HAS BEEN FOUND.  BEFORE USING THIS
  217. C   ROUTINE, THE VECTOR B SHOULD BE SCALED IF MATRIX A WAS SCALED, USING
  218. C   THE SAME SCALE FACTORS.  WITHIN THIS ROUTINE, THE ELEMENTS OF B ARE
  219. C   REARRANGED IN THE SAME WAY THAT THE ROWS OF A WERE INTERCHANGED, USING
  220. C   THE ORDER VECTOR WHICH HOLDS THE ROW ORDERINGS.  THE SOLUTION IS
  221. C   RETURNED IN X.
  222. C
  223. C ---------------------------------------------------------------
  224. C
  225. C   PARAMETERS:
  226. C
  227. C   RLU       THE LU EQUIVALENT OF THE COEFFICIENT MATRIX
  228. C   B         THE VECTOR OF RIGHT HAND SIDES
  229. C   X         SOLUTION VECTOR
  230. C   N         NUMBER OF EQUATIONS
  231. C   NDIM      FIRST DIMENSION OF A IN THE MAIN PROGRAM
  232. C   IORDER    INTEGER ARRAY OF ROW ORDER AS ARRANGED DURING PIVOTING
  233. C
  234. C ---------------------------------------------------------------
  235. C
  236.       DIMENSION RLU(NDIM,N),B(N),X(N)
  237.       DIMENSION IORDER(N)
  238. C
  239. C ---------------------------------------------------------------
  240. C
  241. C   REARRANGE THE ELEMENTS OF THE B VECTOR.  X IS USED TO HOLD THEM.
  242. C
  243.       DO 10 I=1,N
  244.          J=IORDER(I)
  245.          X(I)=B(J)
  246.    10 CONTINUE
  247. C
  248. C ---------------------------------------------------------------
  249. C
  250. C   COMPUTE THE B VECTOR, STORING BACK IN X
  251. C
  252.       X(1)=X(1)/RLU(1,1)
  253.       DO 50 IROW=2,N
  254.          IM1=IROW-1
  255.          SUM=0
  256.          DO 40 JCOL=1,IM1
  257.             SUM=SUM+RLU(IROW,JCOL)*X(JCOL)
  258.    40    CONTINUE
  259.          X(IROW)=(X(IROW)-SUM)/RLU(IROW,IROW)
  260.    50 CONTINUE
  261. C
  262. C ---------------------------------------------------------------
  263. C
  264. C   NOW GET THE SOLUTION VECTOR.  X(N)=B(N) ALREADY.
  265. C
  266.       DO 70 IROW=2,N
  267.          NVBL=N-IROW+1
  268.          SUM=0
  269.          NP1=NVBL+1
  270.          DO 60 JCOL=NP1,N
  271.             SUM=SUM+RLU(NVBL,JCOL)*X(JCOL)
  272.    60    CONTINUE
  273.          X(NVBL)=X(NVBL)-SUM
  274.    70 CONTINUE
  275. C
  276.       RETURN
  277.       END
  278.